Time dependence of critical behavior in multifragmentation 
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We study signatures of critical behavior in microscopic simulations of small, highly excited Lennard- 
Jones drops. We focus our attention on the behavior of the system at the time of fragment formation 
(which takes place in phase space) and compare the results with the corresponding ones obtained at 
asymptotic times (experimentally accessible). The four critical exponents (r,/3, a and 7) found at 
fragmentation time have shown to be stable against time evolution, indicating that the asymptotic 
stage reflects accurately the physics at fragmentation time. Even though evidence of critical behavior 
arises from the calculations, we can not affirm that the system is performing a second order like 
phase transition. 



The process of rapid expansion and fragmentation of 
highly excited drops attracts the attention of physicists 
in different areas, in particular in the field of nuclear 
physics. The possibility of facing a phase transition 
was initially triggered by the work of the Purdue Group 
101 when the resulting fragment mass spectra of proton- 
nucleus colission was fitted by a power law like depen- 
dence. This has initiated a series of works focusing on 
the extraction of signals of critical behavior from both, 
experimental measurements as well as data resulting from 
numerical simulations |^-|^ . In the experimental case the 
only information available corresponds to the asymptotic 
regime and the primordial fragments have to be recon- 
structed from this data by the application of different 
approximations |^ . In a series of works (see for example 
I pof , prjl) critical exponents have been calculated from 
the analysis of the experimental fragment mass distri- 
butions. On the other hand, in the case of numerical 
simulations all the microscopic information is available 
at all times, and if appropriate fragment recognition al- 
gorithms are used, one can unveil the complex fragmen- 
tation mechanism. In this work we use the already pre- 
sented Early Cluster Recognition Algorithm (ECRA) to 
recognize fragments . The advantage of this method- 
ology is that it is able to find, very early in the evolution, 
the partitions of the system that give rise to the asymp- 
totic fragments. These partitions (ECRA clusters) are 
the most bound density fluctuations in phase space and 
has be en pro ved to be the seeds of the asymptotic frag- 
ments [l^jlj]. As all order correlations are available, it 
is possible to follow the evolution in time of the ECRA 
clusters and determine the time of fragment formation 
{Tff) as the time at which these clusters attain micro- 
scopic stability. Once the distribution of ECRA clusters 
at Tff is known, we can look for signals of critical be- 
havior and extract the values of the critical exponents r, 
/9, tJ and 7 (for a definition of these critical exponents 
see section III) . We then test the stability of the values 
of these exponents as a function of time, checking if the 
values extracted at asymptotic times are similar to the 



ones obtained at Tff. In this way, we try to determine 
if the asymptotic stage information properly characterize 
the state of the system at fragmentation time. 

This work is organized as follows: In Section I we de- 
scribe the model used to simulate the fragmentation pro- 
cess as well as the algorithms adopted to recognize the 
fragments. In Section II we will give a brief description of 
the scaling model applied to analyze the data extracted 
from simulations, as well as the behavior of the system 
near the critical point. In Section III we analyze the be- 
havior of some signals of critical behavior at Tf f and at 
asymptotic time. The selected signals of critical behav- 
ior are: the second moment of the fragment distribution 
{M2), the normalized variance of the mass of the max- 
imum fragment {NVM), the best power law fit to the 
fragment spectra according to the assumed scaling model 
(minimum x^) and the critical exponents r, 7, (3 and a. 
Finally, in Section IV, conclusions are drawn. 



I. COMPUTER EXPERIMENTS 



In previous works [|13[-|15[ we have already analyzed the 
dynamics of fragment formation and the caloric curve for 
excited drops made up of 147 particles interacting via a 
6-12 Lennard Jones potential. The interaction potential 
reads: 



V{r) 




(f)=-(ft)'^(ft)' 



r < Tc 
r > r„. 



(1) 



We took the cut-off radius as rc = Scr. Energy and dis- 
tance are measured in units of the potential well (e) and 
the distance at which the potential changes sign (cr), re- 
spectively. The unit of time used is: to = y/ahnJiSe. In 
our numerical experiments initial conditions where con- 
structed using the already presented |jl^,|l^ method of 
cutting spherical drops composed of 147 particles out of 
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equilibrated, periodic, 512 particles per cell L.J. system. 
One of the key ingredients in the analysis of fragmen- 
tation is the determination of the time at which frag- 
ments are formed. This can be accomplished if a proper 
fragment recognition method is employed. In this pa- 
per, clusters are defined according to the ECRA algo- 
rithm in which the fragments are related to the most 
bound density fluctuations in phase space, as was men- 
tioned above. These fragments are referred as ECRA 
clusters. In this case fragments are given by the set of 
clusters {C^} for which the sum of the fragment internal 
energies attains its minimum value: 

^^fn* = [E^r-+ E ^.-.fc] (2) 

jec, j.keCij<k 

where the sum in the first equation of (j^) is over the 
clusters of the partition, Kj "^- is the kinetic energy of 
particle j measured in the center of mass frame of the 
cluster which contains particle j, and Vij stands for the 
inter-particle potential. 

Another way of recognizing fragments is provided by 
the Minimum Spanning Tree (MST) algorithm ||l|,|l§. 
In this case, given two particles they belong to a 
cluster C if the following relation is satisfied: 

ieC^3jeC/nj <r,i (3) 

with r,jj the intcrparticle distance and the cluster- 
ization radius ( rd < rc = 3a). In our case, we took 
fci = fc- It should be noticed that in this definition the 
correlations in momentum space between particles are 
completely disregarded. 

The asymptotic fragments detected by MST algorithm 
are observables. On the other hand, the ECRA clus- 
ters, being defined in phase space, are observables only 
when they coincide with MST clusters. For earlier times 
they are usually embedded into the MST clusters. It 
has been shown ||l^-|l^,|l^ that different properties of 
the ECRA fragments (for instance the average size of 
the maximum fragment, or the average multiplicity of 
fragments in given bins, etc.) become stable very early 
in the evolution, much earlier than the times of stabiliza- 
tion of the same quantities calculated from MST clus- 
ters. As already mentioned in the introduction, the time 
of fragment formation is associated to the time at which 
the ECRA clusters attain microscopic stability. This oc- 
curs when the system switches from a regime dominated 
by fragmentation to one in which the dominant decay 
mode is evaporation. In order to estimate this time (t//) 
we use the so called. Short Time Persistence (STP), in 
which we calculate the stability of the ECRA clusters 
against evaporation-fragmentation and coalescence. It is 
defined in the following way: at a given time t we ana- 
lyze each fragment C* of size by searching on all the 
fragments (7*+"^* present at time t + dt for the biggest 
subset iVmax]*^'** of particles that belonged to C*. Then 



we assign to this fragment a value STP^ = ""^y — . 
This term account for the " evaporation- fragmentation 
process", but on the other hand one has to take care 
of the cases in which the A^max]*^''* does not consti- 
tute a free cluster but is embedded in a bigger frag- 
ment of mass iV*'*''**. We include this effect by defining 

ST Pi ~ '^t+dt ■ Finally the Short Time Persistence 
reads : 

STPit,dt) = ^(^^lMl^^^I^^) ) (4) 

where (...)„j is the mass weighted average over all the 
fragments with size iV > 3. And (...)g is the average over 
an ensemble of fragmentation events at a given energy 
E or multiplicity m. From this definition it is clear that 
STP ~ 1 if the analyzed partition remains unaltered dur- 
ing dt and, on the other hand, STP if its microscopic 
composition presents large changes during that lapse of 
time. 

In Fig. 1 we show one typical calculation of this mag- 
nitude. This is for a very energetic case E = 2.65e for 
which a fast exponentially decaying spectrum is obtained. 
In this figure we show STP{t, dt) as a function of t and for 
different values of dt, namely dt = 0.5to, lio-j l-5^0j 2to. 
The almost horizontal line denotes the reference value 
of STP and corresponds to fragments undergoing only a 
simple evaporation process. In this way we say that when 
STP crosses this line the behavior of the system goes, on 
the average, from fragmentation to evaporation, and we 
call this time Tf f . Using the above mentioned criteria we 
have obtained time-scales ranging from Tff 75to for 
E = -0.5e up to T// - lOio for E = 3.7e. 
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time (t„) 



FIG. 1. Short time persistence (STP) as function of time 
for four different values of dt (empty squares:cit — 0.5to, full 
circles Ito, empty triangles 1.5to and stars 2to) and the refer- 
ence value (solid line) 

In the same way, but working with the fragments given 
by the MST algorithm we get the times of fragment emis- 
sion Tfe. ,i.e. the times at which the clusters defined 
in q-space become stable. These times are much larger 
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than the corresponding time of fragment formation [rff). 
For example, for E = — 0.5e, r/e = lOOio- We took as 
asymptotic time, (tq), a value much larger that all the 
Tfe calculated, i.e., — 600to. 

Once we determine the time of fragment formation 
{Tf f) and the asymptotic time (tq), it is possible to calcu- 
late the corresponding ECRA fragment mass distribution 
at Tf f and compare it with the corresponding MST frag- 
ment mass spectra at Tq . The result of such a calculation 
are displayed in Fig. 2. Three typical mass spectra are 
shown: U-shaped in panel (a), power law like in (b) and 
exponential decaying one in (c). Filled circles denote the 
ECRA distributions at time of fragment formation and 
empty triangles the MST ones at asymptotic times. It 
can be notice that the resulting distributions are esen- 
tially the same. 
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FIG. 2. Three typical fragment distributions {ua' number 
of clusters per particle of size ^4) corresponding to three dif- 
ferent excitation energies and classified by their multiplicity 
values. Black circles correspond to ECRA fragment distribu- 
tions at Tf f and empty triangles to MST ones at asymptotic 
times. The values of m indicate the multiplicity of the cor- 
respondig mass spectra . In panel (a) we show the U-shaped 
distributions (low excitation energies), in (b) the power law 
and in (c) the exponential decaying ones (high excitation en- 
ergies) 



II. THE SCALING MODEL: CRITICAL 
EXPONENTS AND OTHER SIGNATURES OF 
CRITICALITY 

The behavior of the system near its critical point can 
be described in terms of scaling laws. A widely used scal- 
ing model for the fragment mass distributio n ilO| ,[rT|,Pl , 
based in renormalization group arguments |18[, can be 
applied to this kind of problem. In this approach, the 
cluster distribution can be written as: 

UAie) = qoA-^f{z), (5) 

where ha is the probability per particle of having a 



cluster of size A, r is a critical exponent, qo a normaliza- 
tion constant and f{z) a scaling function. This function 
depends on the distance to the critical _2pint, e, and the 
mass number, A, via the combination 



(6) 



with a a critical exponent. The distance to the crit- 
ical point can be defined as e = {Tc — T)/Tc for usual 
thermodynamic systems, and e — [nic ~ m)/nic [|lO| for 
multifragmentation experiments. 

At the critical point, e = (with /(O) = 1 [|l9|), and 
the cluster mass distribution becames a power law: 



(7) 



where the normalization go depends on r via a Rie- 
mann C, function pOf: 



1 



90 = 



T.aA^-^' 



(8) 



This dependence, coming from the normalization con- 
dition Mi(e = 0) — 1, should be taken into account to 
properly extract the critical exponent r when fitting the 
cluster mass distribution to a power law. 

A relevant quantity in the analysis of critical behavior 
is the second moment of the mass distribution: 



(9) 



It has been shown (see for example [|T^) that M2 di- 
verges at the critical point if the fragment distribution 
satisfies equation (^). This divergence can be described 
in terms of the critical exponent 7: 



M2(e)cxr±|e| 



(10) 



Another relevant magnitude in the search of critical 
behavior is the normalized variance of the size of the 
biggest cluster [NVM) pll, defined by: 



NVM = 



< A, 



< A,, 



< A,, 



> 



(11) 



where Amax is the size of the biggest fragment and 
< ... > is the average over an ensemble of fragmentation 
events with a given multiplicity m. In it was shown 
that it is a robust signature of critical fluctuations. 

It is worth noting that for finite systems it is expected 
that scaling assumptions would be valid in a range of 
masses where finite size effects can be avoid and where 
the behavior of the system resembles the behavior of the 
infinite one. 

We will see in the next section that for exploding LJ 
drops, the fragment distributions at Tff and Tq follows 
the scaling model within a rather broad range of mass. 
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III. CRITICALITY ANALYSIS 

In what follows we search for critical behavior in our 
numerical simulations via the analysis of the fragment 
mass disttributions at both, fragmentation time (t^/) 
and asymptotic times (tq). 



FIG. 3. Normalized variance of the size of the maximum 
fragment (NVM) (a), second moment of the cluster distribu- 
tion (M2) (b), and best fit to a power law (with qo = ?o(t) 
via eq.(^)) (c) as function of the multiplicity for fragmenta- 
tion time (Tff), full triangles, and for asymptotic time (ra), 
empty triangles. In (d) we can see average multiplicity and 
its standard desviation for each energy for t// (fuU circles) 
and for Ta (empty circles) 



A. Signals of critical behavior 

We first calculate the normalized variance of the mass 
of the biggest fragment, NVM, and the second moment 
of the cluster mass distributions, M2 as a function of 
the multiplicity, being the multiplicity a relevant observ- 
able in multifragmentation. The biggest fragment of each 
event should be removed for the liquid phase (m < m^) 
in the calculation of M2 Because we don't know a 

priory the value of the critical multiplicity, the biggest 
fragment was excluded for all events . 

The dependence of these magnitudes with the multi- 
plicity is shown in Fig. 3 for the two times analyzed, 
Tff and Ta- It can be seen that both observables peak 
at, essentially, the same multiplicity value for each of the 
analized times. At Tff, the NVM peaks for m = 60 ± 1 
[Fig. 3(a)] and the M2 for m = 65 ± 4 [Fig. 3(b)] (black 
triangles in both cases), whereas at asymptotic times, 
the magnitudes peak for to = 73 ± 1 [Fig. 3(a)] and 
TO = 77±3 [Fig. 3(b)] respectively (empty symbols). The 
higher multiplicity values at asymptotic times are due to 
evaporation processes undergone by the fragments, which 
increase the number of monomers and dimers. If we con- 
sider only fragments larger than A = 2, both magnitudes 
{M2 and NVM) peaks at the same value of multiplic- 
ity at both times. These results are compatible with the 
presence of critical behavior. 
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We will now search for a power law mass spectra ac- 
cording to the scaling model depicted in equation (||). 
As we have seen, the normalization Mi(e = 0) = 1 im- 
poses a condition on qo such that qq = qo{T) via a Rie- 
mann ( function (eq.^. Therefore, if we want to find a 
power law behavior in the mass spectra obtained from 
the numerical simulations, we must fit them with equa- 
tion (^ taking into account that the normalization, qo, 
is given by equation (^). This method was successfully 
used in [|l^ both, for percolation and Au + C multifrag- 
mentation. In our case, the best fitted mass spectra by 
/(A, r) = qo{T)A~'^ determines the critical multiplicity 
whereas the corresponding slope in a log-log plot, the 
critical exponent r. The quality of the fitting procedure 
was measured via the standard coefficient (i.e., a min- 
imum in corresponds to the best fit). 

The fitting was performed in the mass range Q.02Atot < 
A < 0.15 Atot in order to avoid finite size effects. In 
Fig. 3(c) we can see the behavior of the coefficient 
as a function of the multiplicity as calculated at r/ / (full 
triangles) and Ta (open triangles). The minimum in 
signals the critical multiplicity: nic = 59 ± 2 for r/ / and 
TOc = 76 ± 2 for Ta ■ (The obtained values of r from these 
fitting procedures will be reported in next sub-section). 
These results agree very well with those obtained above 
via NVM and M2. 

The relation between the multiplicity and the energy 
of the system is plotted in Fig. 3(d). We can see that 
the critical multiplicities correspond mainly to energies 
around E ~ 0.3e, energy for which the collective motion 
begins to be noticeable [0. It can be seen that there is 
a one to one correspondence between the energy and the 
multiplicity and, moreover, starting right bellow rUc the 
relation between energy and multiplicity becomes almost 
linear. This feature supports the use of multiplicity as a 
control parameter. A similar behavior was also observed 
in percolation for the relation between the multiplicity 
and the bond probability p. 

All these signals show evidences of critical behavior 
not only at asymptotic times, when fragments are ex- 
perimentally accessible, but also at fragmentation time, 
when the particles of the system are still interacting and 
they are recognized by the ECRA algorithm. 
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B. Critical exponents 



rewritten as: 



We now calculate four critical exponents: t, 7, cr, and 
(3 at time of fragment formation (rff) and at asypmtotic 
time (to). 

The calculation of r is performed according to the pro- 
cedure described in the previous section, i.e., finding the 
fragment mass spectra that is best fitted by a power law 
according to the scaling model in a given range of masses. 
From this calculations we got r — 2.18 ± 0.03 for both 
times, Tff and r^. The biggest fragment of each event 
was excluded when generating the average cluster distri- 
bution of a given multiplicity We can see that the 
resulting value of r does not depend on the fit range as 
long as we stay in the range 0. 02 Atot < A < 0.15 Atot, 
where finite size effects are not important. The cluster 
mass distributions at the critical multiplicity and their 
corresponding fits arc shown in Fig. 4(a) for Tff and 
Fig. 5(a) for Ta- 



r.A- 



(13) 




FIG. 4. The critical exponents at r//; r (a), a (b) , (3 (c) 
and 7 (d). In (a) we plot the number of clusters per particle, 
n{A), as a function of the mass number A and its best fit 
(full line). In (b), the value of emax (defined in the text), 
is plotted against the mass number A and its best fit (full 
line). The slope gives the value of a. In (c) we plot the mass 
of the biggest fragment normalized to the size of the system 
(Amax/Atot) as a function of e = and its fit giving the 

value of (3. In (d) we plot the M2 as a function of e and the 
curve with the value of 7 obtained in the fitting procedure 
(full line). 

We have seen in the previous section that the argument 
of the scaling function, z, can be expressed in terms of 
a via equation (^. Taking into account that f{z) has a 
single maximum, we can calculate a in terms of the multi- 
plicity by looking for the value of e (emax) that maximizes 
the production of clusters of a given size A: 



i^max) QqA f {z^dx^ 



(12) 



Therefo re, by plotting emax , as a function of A we get a 
and z max |10|2^. The results obtained are ct = 0.51±0.15 
for Tff and a = 0.64±0.18 for Ta, and they are plotted in 
figures 4(b) and 5(b) respectively. Although the values 
of a at time of fragment formation and at fragmentation 
time do not agree as well as r, they do agree when error 
bars are considered. 

We calculate then the (3 exponent related to behavior 
of the order parameter, which in our case is the biggest 
fragment mass {Amax)- The behavior in the limit e — > 
in an infinite system is given by: 



(^maxMtot) oc e < 0, 



(14) 



where A„iax is the average mass of the biggest fragment 
at a given multiplicity, At^t is the mass of the system and 
e the distance to the critical point as defined above. In 
finite systems, ^ is valid in an intermediate region of 
e where finite size effects are negligible. The relations 
(A-max/Atot) vs 6 are displayed in figures 4(c) for Tff 
and 5(c) for Ta- The obtained values of /3 are 0.29 ± 
0.08 and 0.28 ± 0.13 respectively. These values are very 
close to 0.31, which is the one corresponding to liquid- 
gas transition (3D-Ising Universality class ) , and far from 
the percolation one (0.45). However, we beheve that the 
smallness of the system is responsible of the large errors 
in /? and of a possible sub- valuation of the exponent, as 
was noticed in |22|. 
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time Ta 



where the argument of / is Zmax = A'^ e^nax- This can be 



As we have shown in equation (^^ , the behavior of M2 
near the critical point can be described in terms of the 
critical exponent 7. However the divergence predicted in 
this equation for the limit e — > 0, is valid for an infinite 
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system. For a finite system, we have to find an inter- 
mediate region of e where M2 cc e"' (In this calculation, 
the biggest fragment of each event is removed when the 
multiplicity is below TOc). The procedure used to deter- 
mine 7 is the following first, we take a value of TOc 
from the range given by NVM , M2 and the best fitted 
mass spectrum. Then we determine two ranges of fitting 
boundaries, one for the gas region (e < 0, [e™™ : e™°^]) 
and the other for the liquid region (e > 0, [g™*" : e™"*^]). 
In each one of these regions, we fit the M2 [m) to a power 
law, saving the exponents 7^ and 7;, and their corre- 
sponding Xg and X? for the gas and the liquid regions 
respectively. For each value of rric we make several fitting 
procedures in a wide range of e, choosing those exponents 
(7g and 7;) that fulfil two conditions: first, the values 
should match each other within the error bars given by 
the fitting routine, and second, the of the fits should 
lie in the lowest twenty percent of the distribution. The 
values that satisfy these criteria were histogrammed and 
the average value (< 7 >) and the variance (cr^) were 
obtained from the corresponding distribution . At r// 
we got 7 = 0.77 ± 0.25 and at r^, 7 = 0.72 ± 0.33. The 
results obtained are displayed in figures 4(d) and 5(d) for 
T/ / and Ta respectively. We can see that the value of 7 
is smaller than the one expected for a liquid-gas phase 
transition, although far away from the percolation value 
(see table I). 

In Table 1 we can see the values of the critical expo- 
nents calculated at time of fragment formation (rjy) and 
at asymptotic times {to). Moreover, the values of these 
exponents for two different universality classes, 3D-Ising 
and Percolation, are also presented. It is clear from this 
table that: 

1 - The time evolution after fragmentation time does 
not alter significatively the values of the critical expo- 
nents, i.e. the exponents obtained from asymptotic frag- 
ments distributions properly characterize the fragmenta- 
tion process. 

2 - The exponents are closer to the liquid-gas univer- 
sality class than to the percolation one. However, despite 
the strong evidence arisen in the previous subsection (see 
Fig. 3), we can not be conclusive about the critical be- 
havior of the drop: the error bars in a avoid us to discern 
if the transition could be cast in some of the universality 
clasess presented; the value of (3 is closer to 3D-Ising but 
it could be sub-valuated due to the smallness of the sys- 
tem, and finally 7 shows a value close to 3D-Ising, but 
a little bit lower than the expected value for this univer- 
sality class. 

We think that two main features of the studied sys- 
tem should be taken into account in order to analyze 
the obtained results. First its finite size, that smoothes 
any signal of the possible phase transition and, conse- 
quently, induce quite large error bars for the calculated 
exponents. Second, even a certain degree of equilibration 
can be achieved at t//, non-equilibrium effects could still 
be noticeable (see [?]) and be responsible of the disagree- 



ment between the obtained 7 exponent and the expected 
value from the 3D-Ising Universality class. 



TABLE I. Critical exponents calculated at fragmentation 
and asymptotic time, and the corresponding to the SDIsing 
(liquid-gas) and Percolation universality class [20,21] 



Exponents 






3D-Ising 


Percolation 


r 


2.18 ±0.03 


2.18 ±0.03 


2.21 


2.18 


(7 


0.51 ±0.15 


0.64 ±0.18 


0.64 


0.45 




0.29 ±0.08 


0.28 ±0.13 


0.33 


0.41 


7 


0.77 ±0.25 


0.72 ± 0.33 


1.23 


1.82 
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IV. CONCLUSIONS 

In conclusion, we have calculated different signals of 
critical behavior at fragmentation time as well as at 
asymptotic time. We have found that the dynamical evo- 
lution that follows the fragmentation process does not 
change the system criticallity features, i.e., the critical 
exponents values calculated at time of fragment forma- 
tion. In this way, this work shows that it is possible to get 
the critical exponents that describes the physics involved 
in the fragmentation process from the mass distribution 
obtained experimentally. On the other hand, wc can not 
assure than the process under study can be cast into the 
3D-Ising universality class nor into the percolation one. 
In particular, the 7 exponent is far too low from the cor- 
responding value for 3D-Ising class. This might suggest 
that non-equilibrium effects play a crucial role at time 
of fragment formation. Actually, at this time the system 
is in expansion and the collective motion is noticeable 
for excitation energies near the critical point, as was ob- 
served in the caloric curve of Ref. . 
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